{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "IN_COLAB = 'google.colab' in sys.modules\n",
    "if IN_COLAB:\n",
    "    !git clone https://github.com/cs357/demos-cs357.git\n",
    "    !mv demos-cs357/figures/ .\n",
    "    !mv demos-cs357/additional_files/ ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "from PIL import Image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing the SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1) For a square matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 4\n",
    "n = m\n",
    "A = np.random.randn(m, n)\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using numpy.linalg.svd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U, S, Vt = la.svd(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using eigen-decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now compute the eigenvalues and eigenvectors of $A^TA$ as `eigvals` and `eigvecs`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigvals, eigvecs = la.eig(A.T.dot(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Eigenvalues are real and positive. Coincidence?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigvals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`eigvecs` are orthonormal! Check:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigvecs.T @ eigvecs "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now piece together the SVD:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2) For a non-square square matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 3\n",
    "n = 5\n",
    "A = np.random.randn(m, n)\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can obtain the SVD in the full format using `full_matrices=True` (full_matrices=True is the default value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U, S, Vt = la.svd(A,full_matrices=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(U)\n",
    "print(U.shape)\n",
    "\n",
    "print(Vt)\n",
    "print(Vt.T.shape)\n",
    "\n",
    "print(S)\n",
    "print(S.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the eigen decomposition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or you can use get the reduced form of the SVD:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U, S, Vt = la.svd(A,full_matrices=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('A = ', A.shape)\n",
    "print('U = ', U.shape)\n",
    "print('S = ', S.shape)\n",
    "print('V = ', Vt.T.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Relative cost of matrix factorizations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy.linalg as npla\n",
    "import scipy.linalg as spla\n",
    "from time import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_values = np.logspace(1,3.5,10).astype(np.int32)\n",
    "n_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def matmat(A):\n",
    "    A @ A\n",
    "\n",
    "for name, f in [\n",
    "        (\"lu\", spla.lu_factor),\n",
    "        (\"matmat\", matmat),\n",
    "        (\"svd\", npla.svd)\n",
    "        ]:\n",
    "\n",
    "    times = []\n",
    "    print(\"----->\", name)\n",
    "    \n",
    "    for n in n_values:\n",
    "        A = np.random.randn(n, n)\n",
    "        \n",
    "        start_time = time()\n",
    "        f(A)\n",
    "        delta_time = time() - start_time\n",
    "        times.append(delta_time)\n",
    "        \n",
    "        print(\"%d - %f\" % (n, delta_time))\n",
    "        \n",
    "    plt.loglog(n_values, times, label=name)\n",
    "\n",
    "plt.legend(loc=\"best\")\n",
    "plt.xlabel(\"Matrix size $n$\")\n",
    "plt.ylabel(\"Wall time [s]\");\n",
    "plt.grid();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SVD Applications"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1) Rank of a matrix "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creating matrices for the examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 6\n",
    "n = 4\n",
    "# Creating the orthogonal U and Vt matrices\n",
    "X = np.random.randn(m, m)\n",
    "U, _ = la.qr(X)\n",
    "X = np.random.randn(n, n)\n",
    "Vt, _ = la.qr(X)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2) Low-rank approximations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with Image.open(\"figures/quad.jpg\") as img:\n",
    "    rgb_img = np.array(img)\n",
    "rgb_img.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = np.sum(rgb_img, axis=-1)\n",
    "img.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20,10))\n",
    "plt.imshow(img, cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u, sigma, vt = np.linalg.svd(img)\n",
    "print('A = ', img.shape)\n",
    "print('U = ', u.shape)\n",
    "print('S = ', sigma.shape)\n",
    "print('V.T = ', vt.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(sigma, lw=4)\n",
    "plt.xlabel('singular value index')\n",
    "plt.ylabel('singular values')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.loglog(sigma, lw=4)\n",
    "plt.xlabel('singular value index')\n",
    "plt.ylabel('singular values')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write the matrix $A$ as a linear combination of the outer products of right and left singular vectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another approach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k=50\n",
    "compressed_img = u[:,:k+1] @ np.diag(sigma[:k+1]) @ vt[:k+1,:]\n",
    "plt.figure(figsize=(20,10))\n",
    "plt.imshow(compressed_img, cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the error of the low rank approximation?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "error = la.norm(img - compressed_img,2)\n",
    "print(error)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma[k:k+3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Low rank approximation is a method for image compression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "original_size = img.size\n",
    "compressed_size = u[:,:k].size + sigma[:k].size + vt[:k,:].size\n",
    "print(\"original size: %d\" % original_size)\n",
    "print(\"compressed size: %d\" % compressed_size)\n",
    "print(\"ratio: %f\" % (compressed_size / original_size))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with Image.open(\"figures/Foellinger_Auditorium_2007.jpg\") as img:\n",
    "    rgb_img = np.array(img)\n",
    "foellinger = np.sum(rgb_img, axis=-1)\n",
    "\n",
    "def bestk(A, k):\n",
    "    U,S,V = np.linalg.svd(A, full_matrices=False)\n",
    "    return U[:, :k] @ np.diag(S[:k]) @ V[:k, :]\n",
    "\n",
    "low_10 = bestk(foellinger, 10)\n",
    "low_20 = bestk(foellinger, 20)\n",
    "low_50 = bestk(foellinger, 50)\n",
    "\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.subplot(2, 2, 1)\n",
    "plt.imshow(low_10, cmap='gray')\n",
    "plt.title(\"k = 10\")\n",
    "plt.subplot(2, 2, 2)\n",
    "plt.imshow(low_20, cmap='gray')\n",
    "plt.title(\"k = 20\")\n",
    "plt.subplot(2, 2, 3)\n",
    "plt.imshow(low_50, cmap='gray')\n",
    "plt.title(\"k = 50\")\n",
    "plt.subplot(2, 2, 4)\n",
    "plt.imshow(foellinger, cmap='gray')\n",
    "plt.title(\"k = 480\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3) Pseudo-inverse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Square matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 4\n",
    "n = 4\n",
    "# Creating the orthogonal U and Vt matrices\n",
    "X = np.random.randn(m, m)\n",
    "U, _ = la.qr(X)\n",
    "X = np.random.randn(n, n)\n",
    "Vt, _ = la.qr(X)\n",
    "\n",
    "\n",
    "# Creating the singular values\n",
    "S = np.zeros((m,n))\n",
    "# This creates a full rank matrix\n",
    "r = min(m,n)\n",
    "# This creates a rank deficient matrix\n",
    "r = np.random.randint(1,min(m,n))\n",
    "\n",
    "print(\"the rank of A is = \",r)\n",
    "# Completing the singular value matrix Sigma\n",
    "sig = np.random.randint(1,50,r)\n",
    "sig.sort()\n",
    "sigmas = sig[::-1]\n",
    "for i,s in enumerate(sigmas):\n",
    "    S[i,i] = s\n",
    "\n",
    "A = U@S@Vt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rectangular matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 6\n",
    "n = 4\n",
    "# Creating the orthogonal U and Vt matrices\n",
    "X = np.random.randn(m, m)\n",
    "U, _ = la.qr(X)\n",
    "X = np.random.randn(n, n)\n",
    "Vt, _ = la.qr(X)\n",
    "\n",
    "\n",
    "# Creating the singular values\n",
    "S = np.zeros((m,n))\n",
    "# This creates a rank deficient matrix\n",
    "r = np.random.randint(1,min(m,n))\n",
    "\n",
    "print(\"the rank of A is = \",r)\n",
    "# Completing the singular value matrix Sigma\n",
    "sig = np.random.randint(1,50,r)\n",
    "sig.sort()\n",
    "sigmas = sig[::-1]\n",
    "for i,s in enumerate(sigmas):\n",
    "    S[i,i] = s\n",
    "\n",
    "A = U@S@Vt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4) Matrix Norms and Condition number"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Square and non-singular matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 4\n",
    "n = 4\n",
    "# Creating the orthogonal U and Vt matrices\n",
    "X = np.random.randn(m, m)\n",
    "U, _ = la.qr(X)\n",
    "X = np.random.randn(n, n)\n",
    "Vt, _ = la.qr(X)\n",
    "\n",
    "\n",
    "# Creating the singular values\n",
    "S = np.zeros((m,n))\n",
    "# This creates a full rank matrix\n",
    "r = min(m,n)\n",
    "# This creates a rank deficient matrix\n",
    "## r = np.random.randint(1,min(m,n))\n",
    "\n",
    "print(\"the rank of A is = \",r)\n",
    "# Completing the singular value matrix Sigma\n",
    "sig = np.random.randint(1,50,r)\n",
    "sig.sort()\n",
    "sigmas = sig[::-1]\n",
    "for i,s in enumerate(sigmas):\n",
    "    S[i,i] = s\n",
    "\n",
    "A = U@S@Vt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given the SVD of A..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u,s,vt = la.svd(A)\n",
    "s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... determine the eucledian norm of $A$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... determine the eucledian norm of $A^{-1}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... determine the condition number of $A$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Square and singular matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 4\n",
    "n = 4\n",
    "# Creating the orthogonal U and Vt matrices\n",
    "X = np.random.randn(m, m)\n",
    "U, _ = la.qr(X)\n",
    "X = np.random.randn(n, n)\n",
    "Vt, _ = la.qr(X)\n",
    "\n",
    "\n",
    "# Creating the singular values\n",
    "S = np.zeros((m,n))\n",
    "# This creates a full rank matrix\n",
    "r = min(m,n)\n",
    "# This creates a rank deficient matrix\n",
    "r = np.random.randint(1,min(m,n))\n",
    "\n",
    "print(\"the rank of A is = \",r)\n",
    "# Completing the singular value matrix Sigma\n",
    "sig = np.random.randint(1,50,r)\n",
    "sig.sort()\n",
    "sigmas = sig[::-1]\n",
    "for i,s in enumerate(sigmas):\n",
    "    S[i,i] = s\n",
    "\n",
    "A = U@S@Vt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given the SVD of A..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u,s,vt = la.svd(A)\n",
    "s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... determine the eucledian norm of $A$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... determine the eucledian norm of $A^{+}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... determine the condition number of $A$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
